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C D , Abstract. Adaptive Finite Element Method (adaptivity) is known to be an eflective numerical tool for some ill-posed 

■ problems. The key advantage of the adaptivity is the image improvement with local mesh refinements. A rigorous proof of this 
property is the central part of this paper. In terms of Coefficient Inverse Problems with single measurement data, the authors 

' consider the adaptivity as the second stage of a two-stage numerical procedure. The first stage delivers a good approximation of 

the exact coefficient without an advanced knowledge of a small neighborhood of that coefficient. This is a necessary element for 
the adaptivity to start iterations from. Numerical results for the two-stage procedure are presented for both computationally 

■ simulated and experimental data. 
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(— I , 1. Introduction. For the ffi'st time, the relaxation property for the Adaptive Finite Element Method 

(adaptivity) for a class of non-linear ill-posed problems was proved analytically in [16| . The relaxation 
property ensures that the adaptivity is worth to work with. In short, the relaxation is a rigorously derived 
estimate, which shows that the solution computed on a finer mesh is more accurate than the one computed on 
a coarser mesh. Unlike classical Well-Posed problems, this property is not automatic for Ill-Posed problems: 
because of the instability of the inversion in the latter case. The main results of the current paper is Theorem 
I 5.2 (section 5), where a proof, simpler than the one of [TB], is presented. Prior to [H] the relaxation was 

■ observed numerically, rather than analytically, in a number of publications about Coefficient Inverse Problem 
(CIPs), see, e.g. [ll[6l[Il[8l[9l[10l[ni[Tl[15]. 

■ In most theorems of this paper (although not in all of them) we consider only the nonlinear finite 
dimensional case. The infinitely dimensional case would likely result in imposing the well known source 
representation condition, i.e. assuming that the solution belongs to the range of a certain compact operator. 

', The latter cannot be effectively verified. On the other hand, since we are focused on applications of our 

theory to CIPs, then the work in a finite dimensional space is well justified by the fact that we actually work 
with finite elements, the number of which cannot be too large in any practical computation. 

In our analytical derivations throughout the paper we assume that the noise level S is sufficiently small. 
This is both a common and natural assumption in the theory of Ill-Posed problems, especially in the nonlinear 

■ case. Indeed, in principle one can hope to get an accurate solution only if the noise level is small. However, 
if the noise is large, then only a very special procedure, which is designed for a specific problem of interest, 
might or might not deliver an accurate solution. Those procedures cannot be described in the framework of 
Functional Analysis, since each such procedure highly depends on many specifics of a problem of interest. 
On the other hand, such a procedure took place for the second numerical example of this paper, which is 
for experimental data. And noise was quite large in this case: see comments in the beginning of Section 8.2. 
This confirms a commonly known observation that the theory is usually more pessimistic than numerical 
examples. 

This paper summarizes recent results of the authors on the relaxation property for the adaptivity for 
Ill-Posed problems, see [TTl [111 . First, results are formulated in the Functional Analysis setting. Next, 
they are applied to a CIP for a hyperbolic PDF. Both formulations and proofs of almost all theorems are 
modified here, compared with above publications. It is shown in section 8 (Remark 8.2) that the relaxation 
property helps to work out the stopping criterion for mesh refinements. Theorems 5.3 and 5.4 as well as 
the numerical example of Test 1 were not published before. 

The essence of the adaptivity consists in the minimization of the Tikhonov functional on a sequence 
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of locally refined meshes. It is important that due to local rather than global mesh refinements, the total 
number of finite elements is rather moderate. If this number would be very large, then the corresponding 
space of finite elements would effectively behave as an infinitely dimensional one. However, in the case of a 
moderate number of finite elements, this space effectively behaves as a finite dimensional one. Since all norms 
in finite dimensional spaces are equivalent, then we use the same norm in the Tikhonov regularization term 
as the one in the original space (except of Section 2.1). This is obviously more convenient for both analysis 
and numerical studies than the standard case of a stronger norm in this term [21 [11] [27j \A5\ [46] . Numerical 
results of the current and previous publications confirm the validity of this approach. Note that although 
the finite dimensional version of the original ill-posed problem might be well posed, at least formally, in 
the actuality it inherits the ill-posedness at certain extent. Thus, the use of the regularization term is still 
important for the stabilization. 

Recall that a minimizer of the Tikhonov functional, if it exists, is called regularized solution of the 
corresponding equation [2l [TTl [23l [27l [45l [46]. It is well known, however, that Tikhonov functionals for 
nonlinear Ill-Posed problems, such as, e.g. CIPs, suffer from the phenomenon of multiple local minima and 
ravines. Hence, many regularized solutions might exist. In addition, there is no guarantee that a gradient- 
like or a Newton-like method of minimizing such a functional would converge to the exact solution x* , unless 
the first guess xq would not be sufficiently close to x* . In other words, those are locally convergent methods, 
so as the adaptivity is. Therefore, the assumption in some theorems of this paper that the norm ||.to — x*\\ 
is sufficiently small is a natural one, and the goal of the adaptivity is to refine xq. 

Assuming that the norm ||a;o — a;*|| is sufficiently small, we estimate below the distance between a 
regularized solution and the one obtained after adaptive mesh refinements. Next, we estimate the distance 
between the latter solution and x*. These are the so-called "a posteriori error estimates" (Theorems 5.2, 
5.3, 6.4 and 6.5 below). This is a new element here. Indeed, in the past publications about the adaptivity 
for ill-posed problems, a posteriori error estimates were obtained only for either the Tikhonov functional or 
the Lagrangian, rather than for solutions themselves, see, e.g. [U [U [HI [3 [H [IHl [131 [II]- 

It follows from the above discussion that, prior to applying the adaptivity to a CIP, it is necessary to 
figure out at least one point in a small neighborhood of the correct solution. Hence, we have developed 
a two stage numerical procedure for some CIPs for a hyperbolic PDE. On the first stage, the so-called 
"approximately globally convergent method" [TTl [HI [HI |32l [35l [ST] [38] delivers the key ingredient of any 
locally convergent method: a good approximation xq for the exact solution x* . On the second stage, the 
adaptivity uses this approximation as a starting point for a refinement [TTl 1131 E] [15] . 

The adaptivity for an ill-posed problem, specifically for a CIP for a hyperbolic PDE, was first proposed 
in 2001 in [6]. Also, in 2001 a similar idea was proposed in |4], although an example of a CIP was not 
considered in [4]. In both these first publications the so-called "Galerkin orthogonality principle" was used 
quite essentially. The adaptivity was developed further in a number of publications, where it was applied 
to CIPs [31 [S] [3 [5] [H [13 • A posteriori error estimates in an approximately globally convergent method 
was derived and an adaptive globally convergent method was developed at the first time in [1]. In [36] a 
posteriori error estimates was presented and an adaptive finite element method was applied for the solution 
of a Fredholm integral equation of a first kind. We also refer to [26] where the adaptivity was applied to a 
parameter identification problem. In a CIP an unknown coefficient of a PDE should be reconstructed using 
boundary measurements. In a parameter identification problem an unknown coefficient is reconstructed 
assuming that the solution of the corresponding PDE is given either everywhere inside of the domain of 
interest or on a grid inside of this domain. In the recent publication [40] the adaptivity was applied, for the 
first time, to the classical Cauchy problem for the Laplace equation and quite accurate images were obtained. 
Unlike other works on this topic, both lower and upper error estimates were obtained in [40j . 

In the sections 2-5 we use the apparatus of the Functional Analysis to address above items 1-4 for rather 
general ill-posed problems. In section 6 we deduce from sections 2-5 some results for a CIP for a hyperbolic 
PDE. In section 7 we present mesh refinement recommendations. In section 8 we present numerical results, 
including ones for real experimental data. In numerical studies of this paper we use the above mentioned 
two-stage numerical procedure. 
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2. Minimizing Sequence and a Regularized Solution Versus the First Guess. In this section 
we estimate the distances between terms of the minimizing sequence of the Tikhonov functional and the 
exact solution via the distance between the first guess and the exact solution. In the finite dimensional case 
the minimizing sequence is replaced with the regularized solution. 

2.1. The infinitely dimensional case. Let B,Bi,B2 be three Banach spaces. We denote norms in 
these spaces respectively as |M| , |Mli , IMI2 • As it is conventional in the theory of Ill-Posed problems, we 
assume that Bi C B, \\x\\ < C||a;|j-^ ,Va; S Bi and Bi ~ B, C ~ const. > 0, and the closure Bi is in the 
norm ||-|| . Furthermore, we assume that any bounded set in Bi is a compact set in B. Let G C Si be a 
set and G be its closure in the norm ||-|| . Let F : G B2 he a one-to-one operator, which is continuous in 
terms of norms ||-j| , |1-||2 . Consider the equation 

F{x) =y,x€G. (2.1) 

As it is usually done in the regularization theory [21 [TTl |23l Ull ESJ |46j , we assume that the right hand side 
of equation ()2.1|) is given with a small error S G (0, 1). We also assume that there exists an "ideal" exact 
solution X* of ()2.ip with the "ideal" exact data y* (in principle, there might be several exact solutions). 
Thus, we assume that 

F{x*)=y*,x* eG,\\y-y*\\^<d. (2.2) 

Let xq G _Bi be a first guess for the exact solution x* . Usually one assumes that xq is located in a small 
neighborhood of x* . Consider the Tikhonov functional 

M„ (x) = l\\F{x)-y\\l + ^ \\x - xoWl - 2:, G G, (2.3) 

where a G (0, 1) is the regularization parameter. We impose a rather conventional assumption that 

a^a{S)^ S^'',^i = const. G (0,1/2). (2.4) 
The second term in the right hand side of (|2.3p is called "the Tikhonov regularization term" . Let 

ma = infM„(x). (2.5) 



Hence, there exists a minimizing sequence {x^j^^ C G such that lim„_^oo Ma {x") = nia. By (|2.2|) . ()2.3p 
and (|2?5|) 

rria < Ma (x*) < 6^ +a\\xo - x*\\l . (2.6) 

Hence, there exists an integer N — N {5,F) > 1 such that Ma (^Xr}^^^ < S'^ + a \\xq — x*\\-^ ,\fn> N. Hence, 
by (EJ) 



1 /2 

^<\/2(<52(i-'') + |l.To-x*||?) +\\xo\\,,yn>NiS,F). (2.7) 



Suppose that an a priori upper estimate of the distance between the first guess and the exact solution is 
given, 

\\xq - x*\\^ < A,A = const. > 0, (2.8) 
"(-5) 



where the number A is given. Then ()2.7p implies that 
P {xq , A) defined as 



< V2{A + l) + \\xo\\^ . Consider the set 
P{xo,A) ^ {x e G : \\x\\^ < V2(A+ 1) + llxolli} • (2-9) 



Let P := P {xq, A) be its closure in terms of the norm || • || . Hence, P C G. Since the set P {xq, A) is bounded 
in terms of the norm ||-j|j^ , then P is a closed compact set in the space B. Consider the range F (P) C P2 of 
the operator F on the set P. Since the operator F : G B2 is continuous in terms of norms ||-|| , ||-||2, then 
F (P) is a closed compact set in B2- Furthermore, since F is one-to-one, then by the foundational theorem 
of Tikhonov [TTl [27l |45l ES] the inverse operator F^^ : F (P) P is continuous. Therefore, there exists the 
modulus of the continuity of the operator F~^ on the set F (P) . This means that there exists a function 
ujp (z) , z S (0, 00) such that 



(z) > 0, ujp (zi) < LOp (^2) if zi < Z2, lim ujp (z) = 0, 

2->0+ 

\\X1 -X2\\< LJp (IIP (Xl) - F {x2)\\2)yXl,X2 £ P. 



(2.10) 
(2.11) 



Theorem 2.1 compares the distance ||.to — 2:*||-^ with the distance between terms of the minimizing 
sequence and the exact solution x*. 

Theorem 2.1 (rate of convergence). Let B,Bi,B2 be Banach spaces, G C Pi be a convex open set and 
P : G — > P2 be a one-to-one continuous operator in terms of norms \\-\\ ,||-||2- Let conditions i2.2\) . (2^, 
\2. 5\] and 112. 8\} be in place. Then for any number S G (0, 1) there exists an integer N ~ N (S,F) > I such 
that 



2.0(5) <ojp(2S''VATtj ,Vn>N{S,F). 



(2.12) 



Proof. Using ([2^ and we obtain for n> N{S, F) 

- P ix*)\l = \\f -y + y-F (x*) 



< 



Fx. 



MS) 



+ \\y-y*\\2 < 



2M„ 



1/2 
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(2.13) 



< ((5^ + S^^' \\xo - x*\\iy^^ + (5 < 2(5'^ (1 + \\xo - a;*||i)^^^ < 2S''y/A + T 



By (EU, dm) and (EH) x* e P. Therefore, ([2ll|) and ((2l3| imply ([2T2)) . □ 



via the distance ||a;o 



between the first guess 



Theorem 2.1 estimates the distance 

and the exact solution for any S G (0, 1) . Still, it is natural to ensure that the distance between terms of 
the minimizing sequence and the exact solution is strictly less than the distance between the first guess xq 
and the exact solution. This can be ensured only for sufficiently small values of the noise level (5. Although 
Corollary 2.1 has a similarity with the well known convergence theorem of the minimizing sequence for the 
Tikhonov functional (see, e.g. page 33 in |11|). still in that theorem only a subsequence converges rather 
than the entire sequence. Besides, estimate (|2.14p is useful by its own right, and also the convergence rate 
(|2.12p . from which (|2.14p is derived, seems to be new. 

Corollary 2.1. Let B,Bi,B2 be Banach spaces, G C Pi be a convex open set and P : G — > P2 be a 
one-to-one continuous operator in terms of norms \\-\\ ,||-||2- Let conditions 112. 2\) . ^2.4^ , h2. 5\] and h2. 8]) be 
in place. Let ^ G (0, 1) be an arbitrary number. Assume first that xq ^ x* . Then there exists a sufficiently 
small number (5o = (5o {F^ l^iCt G (0, 1) such that 



X 



<e||a:o-x*||,VJe (0,,5o),7i>iV(<5,P). 



In the case xq = x* {2.14^ should be replaced with 

<(,y6e{0,So),n>N{d,F). 



(2.14) 



(2.15) 



In particular, if S = 0, then Sq should be replaced with a sufficiently small number ao G (0, 1) and "S G 
(0,(5o)" should be replaced with a G (0,ao) . 
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Proof. First, let xq ^ x* . By (|2.10l) there exists a sufBciently small number (-P^i^iMiC) ^ (Oi 1) such 
that ujp (2(5'" V'T+T) < i llxo - a;*|| , V(5 G (0, b^) ■ Combining this with ([^1^ . we obtain dUHl)- 

Let now xg = x*. Then again there exists a sufficiently small number (5o (i^, A,/^,^) e (0, 1) such that 
ujp (2(5^ + 1) < Combining this with (pi^ . we obtain □ 

2.2. The finite dimensional case. Consider now the finite dimensional real valued Hilbert space. 
Compared with subsection 2.1, the main new point here is that the minimizing sequence is replaced with a 
minimizer, which exists. This case is of our main interest in the current paper because standard piecewise 
linear finite elements form a finite dimensional space. Unlike the above, we now use the same norm in the 
regularization term as in the original space. This is because all norms are equivalent in a finite dimensional 
space. Nevertheless, since the finite dimensional version of the original ill-posed problem "inherits" the 
ill-posedness, at certain extent, it is still important to use the regularization term for the stabilization. 

Let H and Hi be two real valued Hilbert spaces and dimTJ < oo. Norms and scalar products in these 
spaces denote respectively as |M| , (, ) , IMI2 j (: )2 ■ ^^'^ G C be an open bounded set and F : G ^ H2 he 
a continuous operator. We again consider equations (|2.ip . (|2.2p . where x* £ G,y,y* G Hi- The functional 
Ma (x) in p.3p is now replaced with the functional Jq, (x) , 

J„ {x) = i 11^^ (x) - y\\l + f Ik - -Toll' , x e G,xo e G. (2.16) 
The following lemma follows immediately from Weierstrass theorem. 

Lemma 2.1. Let F be the operator defined above in this section. Then there exists a regularized solution 
Xa € G, 

inf Jq (x) = rnin Jq {x) — Ja {xa) ■ (2-17) 
G G 

Although a similar result is valid for the case when the set G is unbounded, we do not formulate it here 
since we do not need it. The following theorem follows immediately from Theorem 2.1 and Corollary 2.1. 

Theorem 2.2. Let Hilbert spaces H,H2, the set G <Z H and the operator F : G ^ H2 be a one-to- 
one continuous operator. Let conditions i2.S\) . i2.8\) . 112. 16\) and {2.1'T^ be in place. Then for any number 

<5e (0,1) 

\\xa{S) -x*\\< ujF (25''y/ATtj . 

Let ^ e (0, 1) be an arbitrary constant. Then there exists a sufficiently small number 60 = 60 {F,A,n,^) e 
(0, 1) such that for all S G (0,(5o) 

C ||xo - x*|j , if xo 7^ X*, 
^, if xo = X*. 

Ln particular, if S — 0, then 5q should be replaced with a sufficiently small number ao G (0,1) and "S G 
(0, (5o)" should be replaced with a e (0,q;o) • 

3. The Local Strong Convexity of the Tikhonov Functional ()2.16p . In [43| the local strong 
convexity of the Tikhonov functional was established for the ease when the underlying operator F has the 
second continuous Frcchet derivative and the source representation condition is in place. In this section we 
prove the local strong convexity of the Tikhonov functional (|2.16p for the case when the operator F has the 
first continuous Frcchet derivative and the source representation condition is not imposed. 

Let H and H2 be two real valued Hilbert spaces. Let scalar products and norms in them be respectively 
(,),||-|| and (,)2JMl2- Let £{H,H2) be the the space of all bounded linear operators mapping H into 
i?2 and let \\-\\^ be the norm in C {H,H2) . Although we do not assume here that H is finite dimensional, 
we still use the same norm j | x — xo j | in the regularization term in (|2.16p as the one in the original space 
H, rather than a stronger norm as in (|2.3|) . This is again because our true goal is to work in a finite 
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dimensional space of finite elements in the adaptivity (section 1). For any a > and for any x Cz H denote 
Va (x) ~ {z Cz H : \\x ~ z\\ < a} . First, we formulate the following well known theorem. 

Theorem 3.1. |41] . Let G Q H be a convex open set and L : G ^ ^ be a functional. Suppose that this 
functional has the Frechet derivative L' (x) € C {H, M) for every point x € G. Then the strong convexity of 
L on the set G with the strong convexity constant p > is equivalent with the following condition 

{L' (x) -L' {z),x- z)>2p\\x~ zf , Vx, z e G. (3.1) 



Theorem 3.2. Let G Q H be a convex open set and F : G ^ H2 be an operator. Let x* ^ G be an 
exact solution of equation \2.1]) with the exact data y* . Let Vi {x*) C G and let 112. S^) holds. Assume that for 
every x € Vi (a;*) the operator F has the Frechet derivative F' [x) G C [H, H2) . Suppose that this derivative 
is uniformly bounded and Lipschitz continuous in Vi {x*), i.e. 

\\F' {x)\\^<Ni,'ixeVi{x*), (3.2) 
\\F' {x) - F' {z)\\^ < N2 \\x - z\\ ,Wx,ze Vi (x*) , (3.3) 

where Ni,N2 = const. > 0. Let 

a = a{5)^S^'', V(5e(0,l), (3.4) 
/i = const, e ^0, ^ . (3.5) 

Then there exists a sufficiently small number Sq — Sq {Ni., N2, fi) G (0,1) such that for all 6 G (0, (5o) 
the functional Ja{S) i^) strongly convex in the neighborhood Vga^ (x*) of x* with the strong convexity 
constant a/4. In the noiseless case with S = one should replace "5q = Sq (Ni, N2, fJ.) G (0,1)" with 
ao = ccq (-/Vi, iV2) € (0, 1) to be sufficiently small and require that a G (0, ao) ■ 

We refer to [TTJ [TB] for the proof of Theorem 3.2 since it is space consuming. Consider now the finite 
dimensional case. 

Theorem 3.3. Let dimH < 00, G d H be an open bounded convex set, and the rest of conditions of 
Theorem 3.2 holds. Let in h2. 1 6]) the first guess xq for the exact solution x* be so accurate that 

||xo-.T*|l<^. (3.6) 

Then there exists a sufficiently small number Sq ~ 6q(Ni,N2,p) G (0,1) such that for every S G (0, (Jq) 
and for a = a (6) satisfying there exists unique regularized solution Xa(s) of equation \2.1]) on the set 

G. Furthermore, x^(^s) € ^5^^ (a;*) . In addition, the gradient method of the minimization of the functional 
Ja{S) (x) , which starts at xq, converges to Xa^s)- Also, if the operator F is one-to-one on Vi{x*), then 
Xa(S) S ^53f'/3 (2;*) ■ In the noiseless case with 6 = one should replace "Sq = Sq {Ni, N2, p) G (0, 1)" with 
Q^o — Q^o (-^1: -^2) G (0, 1) to be sufficiently small and require that a G (0, ao) . 

Proof. By Lemma 2.1 there exists a minimizer x^^^-j G G of the functional Ja{S)- We have Ja{S) {xa{S)) < 
Ja{S) (2;*) • Also, ~ 2:o|| > ||a;Q(i) ~ 2^*11 ~ ll^^o — . Hencc, using (|2.2p . (|2.16p and (|3.6p . we obtain 

that there exists a sufficiently small number Sq = Sq {Ni,N2, p) G (0, 1) such that for every 5 G (0, 5q) 

- x*|| < + 2 ||x* - .xoll < 5'-^ + = {l + h'-'^^ < ■ I = S'^. 

Hence, Xq,(5) G V^a,, (x*) . Since by Theorem 3.2 the functional J„ is strongly convex on the set (x*) 
and the minimizer Xa(s) G V^s^ (x*) , then this minimizer is unique. Furthermore, since by p.6p the point 
xo G VgSf, (x*) , then it is well known that the gradient method with its starting point at xo converges to 
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Let now the operator F be one-to-one. Let ^ G (0, 1) be an arbitrary number and Xq ^ x* . By Theorem 
2.2 we can choose a smaUer number Sq = Sq {Ni,N2, ^) such that 

W^aiS) < C ll^^o - ,yS e (0,(5o) . 

Hence, p.6p imphes that x^^s) G ^^3^/3 (a;*) . If xq = a;*, then by Theorem 2.2 ^ < ^. Choosing 

$ e (0, , we again obtain that Xa(^s) & ^53^/3 i^*) ■ The noiseless case is similar. □ 

4. The Space of Finite Elements. To prove the relaxation property of the adaptivity, we need to 
introduce the space of finite elements. Let ft C M", n = 2, 3 be a bounded domain. Consider a discretization 
of by an unstructured mesh T using non-overlapping tetrahedral elements in M.^ and triangles in such 
that T ~ Ki, ...^Ki, where I is the number of elements in il, and 

D = UxerK ^ KiU K2... U Ki. 

We obtain a polygonal domain D and assume for brevity that D = ^l. We associate with the triangulation 
T the mesh function h = h{x) which is a piecewise-constant function such that 

h{x) = hK VA' e T, 

where Hk is the diameter of K which we define as the longest side of K. Following section 76.4 of [25] , consider 
piecewise linear functions {cj {x,T)}j^^ C C (17), which are called test functions. Functions {cj (x,T)}^^^ 
are linearly independent in ft. Here, N is the global number of nodes in the mesh T. Let {Ni} be the set of 
nodal points of triangle/tetrahedra K for all K G T. Then 

We introduce the finite element space Vh as 

Vh ^ {v{x) eV -.ve C{n), v\k e Pi{K) \/K e r}, 

where Pi (K) denotes the set of piecewise-linear functions on K with 

V={v{x) -.vix) eL2in)}. 

The finite dimensional finite element space Vh is constructed such that Vh G V . 

Let r be the radius of the maximal circle/sphere inscribed in K. We impose the shape regularity assump- 
tion for all triangles/tetrahedra uniformly for all possible triangulations T which we consider. Specifically, 
we assume that 



(4.1) 



Q-i ^ hx ^ ra2, «!, ^2 = const. > 0, \fK £ T, V T, (4-2) 

where numbers ai,a2 are independent on the triangulation T. Let /imax (2^) and /imin (T) be respectively 
the maximal and minimal diameters of triangles/tetrahedra of the triangulation T. We assume evrywhere 
below that 

''■"'"^^^ <CT,VT (4.3) 



^max (-^) 

for a certain positive constant ct- Obviously, the number of all possible triangulations satisfying ()4.2|) . ()4.3 
is finite. Thus, we introduce the following finite dimensional linear space H, 

H^\JVh{T), VT satisfying 



Hence, 



dimiJ < CX3, He (C (n) n (il)) , d^Je (O) , \ff e H. (4.4) 

In ()4.4p "c" means the inclusion of sets. We equip H with the same inner product as the one in L2 (fi) . 
Denote (,) and ||-|| the inner product and the norm in H respectively, ||/||^ := •= ll/ll : ^/ € ff. 

Everywhere below H is this space. We view the space H as an "ideal" space of very fine finite elements, 
which cannot be reached in practical computations. At the same time, all other spaces of finite elements we 
work with below are subspaccs of H. In particular, this means that we assume without further mentioning 
that (j4.2|) and (j4.3|) are valid for all meshes considered below. 

Keeping in mind the mesh refinement process in the adaptivity, we now explain how do we construct 
triangulations {Tn} as well as corresponding subspaces {M„} of the space H which correspond to mesh 
refinements. Consider the first triangulation Ti with rather coarse mesh. We set Mi := Vh (Ti) C H. Suppose 
that the pair (T„, M„) is constructed after n mesh refinements and that the basis functions in the space M„ 
are {cj {x,Tn)}j^-^ ■ We now want to refine the mesh again. We define the pair (r„+i, M„_|_i) as follows. We 
refine the mesh in the standard manner as it is usually done when working with triangular /tetrahedron finite 
elements. When doing so, we keep (|4.2p . Hence, we obtain both the triangulation T„+i and the corresponding 
test functions {ej (cc, T„+i)}^j^^. It is well known that test functions {ej {x,Tn)}^^^ are linearly dependent 
from new test functions {ej (a;, r„+i)}^+^ . Thus, we define the subspace M„+i as 

M„+i := Span (^{e-j (x, T„+i)}^+'^ . 

Therefore, we have obtained a finite set of linear subspaces {M„}^^-^ of the space H. Each subspace Af„ 
corresponds to the mesh refinement number n, M„+i\M„ 7^ and 

M„ C AI„+i C iJ, n e [I, iV - I] . 

Let / be the identity operator on H. For any subspace M C H, let Pm : H M be the orthogonal 
projection operator of the space H onto its subspace M . Denote for brevity P„ := Pm„ ■ Let /i„ be the 
maximal grid step size of T„. Hence, /i„+i ^ Let be the standard interpolant of the function f G H 
on triangles/tetrahedra of T„, see section 76.4 of [25]. It can be easily derived from formula (76.3) of [25] 
that 

||/-/,!|| <A'||V/||i^(j,)/i„,V/ei/, (4.5) 

where K = if (f2, r, ai, 02) = const. > 0. Since /,( e H,yf e H, then by one of well known properties of 
orthogonal projection operators, 

^ ||/-/,!||, V/eif. (4.6) 

Hence, (|4.5p and (|4.6p imply that with a different constant K ~ K (il, r, ai, 02) > 

1|/-P„/K A'||V/||^^(f,)/i„,V/eif. (4.7) 

Since iJ is a finite dimensional space in which all norms are equivalent, it is convenient for us to rewrite 
(|4.7p with a different constant K = K (fl, r, , ai, 02) > as 

- Pnx\\ s$ K \\x\\ K, V.t e H. (4.8) 

5. RelzLxation. Since we sequentially minimize the Tikhonov functional on subspaces {Af„}^^-^ in the 
adaptivity procedure, then we need to establish first the existence of a minimizer on each of these subspaces. 
In this section the set G and the operator F are the same as in Theorem 3.3, and the functional Ja {x) is the 
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same as in (|2.16p . Theorem 5.1 ensures both existence and uniqueness of the minimizer of the functional Jq, 
on each subspace of the space H, as long as the maximal grid step size of finite elements, which are involved 
in that subspace, is sufficiently small. 

Theorem 5.1. Let conditions of Theorem 3.3 hold. In particular, let the operator F : G ^ H2 be 
one-to-one. Let M Q H be a subspace of H and let V^a^ (x*) H M ^ 0. Assume that \\x*\\ < B, where the 
number B > is known in advance. Suppose that the maximal grid step size h of finite elements of M be 
so small that 

where K is the constant in Furthermore, assume that the first guess Xq for the exact solution x* 

in the functional Ja(S) is so accurate that 113. 6\) is in place. Then there exists a sufficiently small number 
^0 = <^o (-^ij -^2, G (0, 1) such that for every S € (0, Sq) there exists unique minimizer XM.a(S) (z GH M of 
the functional Ja on the set GHM. Furthremore, X]^j ,^(^s) € ^s^i^ (x*)nM. In addition, the functional Ja (x) is 
strongly convex on the set V^a^ (x*) D M with the strong convexity constant a (S) /4. Let Xa(s) S ^53^/3 (x*) 
be the regularized solution of equation \2. 1]) . which is guaranteed by Theorem 3.3. Then the following a 
posteriori error estimate holds 



2 II 

\xM,a(S) - ^a(5)\\ < ^ \\J'a {xM,a(5)) 



Note that since in Theorem 5.1 Vl (x*) C G and V^a,, (x*) n Af 7^ 0, then G n A/ 7^ 0. We do not prove 
this theorem here and refer instead to Theorem 4.9.2 of [TT]; also see Theorem 3.2 of [TB] for a similar result. 

Theorem 5.2 (relaxation). Let Af„ G H be the subspace obtained after n mesh refinements, as described 
in section 4- Let hn be the maximal grid step size of the subspace Mn. Suppose that all conditions of Theorem 
5.1 hold with the only exception that the subspace M is replaced with Mn and the inequality I15.1\) is replaced 
with 

^" ^ 3BN^- ^'-'^ 

Let 6 e (0,(5o) , where the number Sq G (0, 1) is defined in Theorem 5.1. Also, let Vg3^ [x*) H Mi 7^ 0. Let 
Xn € V^sm [x*)C]Mn be the unique minimizer of the functional Ja {x) in h2. 1 6\) on the set GHMn (Theorem 
5.1). Let Xa(S) G 1^53f'/3 (2;*) be the unique regularized solution (Theorem 3.3). Assume that 

Xn ^ Xa(S), (5.3) 

i.e. Xa(s) ^ Mn , meaning that the regularized solution is not yet reached after n mesh refinements. Let 
r] G (0, 1). Then one can choose the maximal grid size hn+i = hn+i {Ni, N2, 6, B, K, rj) € (0, hn] of the mesh 
refinement number (n + 1) so small that 

\\xn+l - Xa{S) II < V \\xn " Xa{S) \\ , (5.4) 

where Xn+i S V^s^ (a;*)nA/„+i is the unique minimizer of the functional \2.16\) on the set Gr\Mn+i. Hence, 

||x„+i - Xa(S) II < rT \\xi - Xa(8) \\ ■ (5.5) 

Proof. In this proof we denote for brevity a{S) :— a. Since Vgs^ (x*) n Mi 7^ 0,A/i C Af„ and 
Vs3^ (x*) C Vi (x*) C G, then {Ygs^ (x*) n M„) C {Vi (x*) n M„+i) ^ 0. Since by Theorem 3.2 the functional 
(|2.16p is strongly convex on the set 1/53^ (x*) with the strong convexity constant a/4, then Theorem 3.1 implies 
that 

^ ||x„+i - a;„(5)||^ < {J'a (Xn+l) - J'a {Xa(5)) ,Xn+l - Xa{S)) . (5.6) 
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Since x^+i is the minimizer on G D Af„+i and Xa is the minimizcr on the set G, then 

{J'a (Xri+l) ,z)=0, Vz e Mn+i; J'a {Xa{S)) = 0. (5.7) 

Relations (|5.7p justify the appUcation of the Galerkin orthogonahty principle [Hin]. By (|5.7p 

{j'a (Xn+l) - J'a {xa{S)) , Xn+1 - Pn+lXa{S)) = 0. (5.8) 

Next, Xn+i - Xa(s) = {xn+1 ~ Pn+iXa(s)) + {Pn+iXa{S} " Xa(s)) ■ Hcncc, (|?^ and ((^ imply that 

^ ||a;„+i - a;Q(5)||^ < (j^ (a;„+i) - (a;Q(5)) ,P«+ia;Q(5) - . (5.9) 
It follows from p.3p that conditions (|5.2p and (|5.3p imply that 

(a;a(5)) II < ^3 II (5.10) 
with a constant TVg = N3 (iVi, ^2) > 0. Also, by (g^ 

||a^a(5) - ^'ri+ia:Q(5)|| < if 1 1 a;„(5) 1 1 /i„+i . (5-11) 
Using the Cauchy-Schwarz inequality as well as p.4p . (|5.10l) and (|5.1ip . we obtain from (|5.9p 

Since by one of conditions of Theorem 5.1 we have an a priori known upper estimate ||a;*|| < B, we now 
can estimate the norm ||2;q;(5)||- Since by Theorem 3.3 Xa(5) S Vgai^ {x*) , then 



a^a(5)|| < \\xa{S) - a;*ll + lla;*|| < -:— + B. 
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Hence, (|5.12p becomes 



\Xn+l - Xa(S)\\ < ( ~ + ^ I ^"+1- (5.13) 



2KN3 f5^^ 



Let 77„ S (0, 1) be an arbitrary number. Since ||a;„ — Xa{s) || 7^ 0, then we can choose hn+i = hn+i S, A, K) g 
(0, hn] so small that 

2KN-, /S^'' 



j2m y—+BjK+i<r]\\xn-Xa{S}\\- (5-14) 

Comparing (|5.14l) with (|5.13p . we obtain the target estimate (|5.4p . □ 

Theorem 5.2 provides an estimate of the distance between points Xn+i obtained via adaptive mesh 
refinement and the regularized solution. We now estimate how far are these points from the exact solution 
x*. Theorem 5.3 follows immediately from Theorem 2.2 and ()5.5p . 

Theorem 5.3. Let conditions of Theorem 5.2 hold. Let 5 e (0,(5o), where the number 5q G (0,1) is 
defined in Theorem 5.1. Then there exists a decreasing sequence of maximal grid step sizes {hkYlXi swc/i 
that 

-a;*|| < i)^ \\xi - Xa{&)\\ + top {25^'^/A + l^ ,fc = l,...,n, (5.15) 

where the number A is defined in h2.l^] and the function up is defined in V2.10\) . 112. In particular, let 
^ € (0, 1) be an arbitrary number. Then there exists a sufficiently small number Si ~ Si (A^i, A'^2, A*, C) S (0, Sq] 
and a decreasing sequence of maximal grid step sizes {/ifc})!^]^ such that for all S G (0,(5i) and for k G [i-,n] 

II *ll ^ fell \\ , f ^\\xQ - X*\\ , a Xq ^ X*, , , 

\\xk+l - X W < T] \\xi - Xa(S)\\ + < ^, ifa;o=X*. ^ ' 
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Since /in+i is the maximal grid step size in the entire domain 17, it seems to be at the first glance 
that Theorems 5.2, 5.3 are about mesh refinements in the entire domain rather than about local mesh 
refinements in subdomains, as it is the case in the adaptivity. Assuming that conditions of Theorem 5.2 
hold, we now show that local mesh refinements are also covered by this theorem. Suppose that the domain 
is split in two subdomains, fi = f^i U il2, ^^i H ~ 0. Assume that the function xq is changing slowly in 
rii and has some "bumps" in fl2. These bumps correspond to small inclusions. It is these inclusions rather 
than slowly changing functions, which are of the main applied interest in imaging. Indeed, those small 
abnormalities model, e.g. land mines, tumors, etc. Hence, it is reasonable to assume that x* is also changing 
slowly in rii. Next, because of Theorem 2.2 and because all norms in H are equivalent, it is reasonable to 
assume that the regularized solution Xa is also changing slowly in fJi.Thus, inequality (j5.17p of Theorem 5.4 
is a reasonable one. Furthermore, it is reasonable to assume that mesh refinements do not take place in 
but only in f22. 

Theorem 5.4 (relaxation for local mesh refinements). Assume that conditions of Theorem 5.2 hold. Let 
h^^^ be the maximal grid step size in il^.Then there exists a sufficiently small number So = So (A^i, A^2,m) ^ 

(0, 1) and a decreasing sequence of maximal grid step sizes |/ifej- such that if the norm \\\7xa{S)\\ji^ ^ 
is so small that with the constant = N3 {Ni,N2) > from I15.10\) 

II lli^(oi) ^ f Ikfc - || , fc = 1, n, (5.17) 

then i5.15\) and Ii5.16\) hold with the replacement of {hkYj^^i w*^'* ^ /^fc f 
Proof. By and ([SAO]) 

||a;fc+l - Xa{S) II < \\^a(S) - Pk+lXa{S) \\ = (5.18) 

2N3 /|| II II II \ 

[\\Xa{S) - Pk+lXa(5) II ^^(f^^) + \\Xa{5} - Pk+lXa{S) \\ ^^(t^^) J ■ 

By gJl) and (|5Tfl) 

2N3 II „ II 2A'iV3 II II , 77 II II , , 

\\x„^s) - Pk+ix^(^5)\\L,(n,) ^ lr^"WllL„(ni)'' - 2 lr^"^«wll- ^^-^^^ 

Next, we obtain similarly with (|5.14p 

II D II ^ II II fK on^ 

\\x^^s) - Pk+iXo.is)\\i^^^^^-^ < ^\\xk-Xo,^s)\\- (5.20) 

It follows from ((5l^ - ((Og|) that 

||a;fc+i -XQ,(i-)|| < 77||xfc -a;„(5)|| ,k = l,...,n. 

Hence, ([53]) holds. Finally, (fSTe]) follows from (|5?5|) and Theorem 2.2. □ 

Remark 5.1. Theorems 5.3 and 5.4 claim that the accuracy of the solution improves with mesh 
refinements, i.e., the relaxation takes place. Comparison of (|5.3p with ()5.15p and ()5.16p shows that the 
solution is adaptively refined until reaching the regularized solution Xa(s)- It is important that by Theorem 
2.2 the accuracy of Xa(s) is better than the accuracy of the first guess xq. Indeed, this ensures that it is 
worthy to work with the adaptivity in order to improve the accuracy of the regularized solution via mesh 
refinements. 

6. Adaptivity for a Coefficient Inverse Problem. We now reformulate some of above theorems 
for the case of a specific GIF. To save space, we do not prove theorems of this section. Instead, we point to 
those results of Chapter 4 of [TT] from which these theorems can be easily derived. 
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6.1. Coefficient Inverse Problem and Tikhonov functional. Let $7 C M"^ be a convex bounded 
domain with tlic boundary dfl S C'^. Let the point xq ^ Q. For T > denote Qx ~ ^7 x {0,T),St = 
dfl X (0,T) . Let d > 1 be a certain number, lu G (0, 1) be a sufficiently small number, and the function 
c (x) G C (R^) be such that 

c (x) E {1 — LL!,d + uj) in fi, c (x) — 1 outside of il. (6-1) 

Below we specify c (x) more. Consider the solution u (x, t) of the following Cauchy problem 

c(x)mm = A^t,a;GM^^e (0,r), (6.2) 
u (x, 0) = 0, ut (.X, 0) = (5 (x — xq) . (6.3) 

Equation (j6.2p governs propagation of acoustic waves, in which case c(x) = 1/6^ (x), where h(x) is the 
sound speed and u{x,t) is the amplitude of the acoustic wave [22]. In addition, (|6.2p governs propagation 
of the electromagnetic field in 2-d, in which case c (x) = Sr (x) is the spatially distributed dielectric constant 
and u (x,t) is one of components of the electric field [44]. Although in the latter application equation (|6.2|) 
is valid only in 2-d, we have successfully used this equation to work with experimental data, which are 
obviously in 3-d, see [TTJ [131121] and section 8. This was explained in Test 4 of [T5]. It was shown in this test 
that the component of the electric field, which was initially sent in a rather simple medium, dominates two 
other components. It was also shown that the propagation of the dominated component is well governed by 
equation (|6.2p . 

Remark 6.1. An alternative to the point source in (j6.3|) is the incident plane wave in the case when 
it is initialized at the plane {xa = xa^o} such that {x3 = X3.0} H f2 = 0. The formalism of derivations below 
is similar in this case. In our derivations below we focus on (|6.3p . because this is the most convenient case 
for derivations. However, in numerical studies we use the incident plane wave, because this case has shown 
a better performance than the point source. 

Coefficient Inverse Problem (CIP). Let conditions (|6.ip - (|6.3p hold. Assume that the coefficient c (x) 
is unknown inside the domain D,. Determine this coefficient for .x G 51, assuming that the following function 
g {x,t) is known 

u\sT=9{x,t)- (6-4) 

The function g (x, t) can be interpreted as the result of measurements of the wave field u (x, t) at the 
boundary of the domain of interest i}. Since the function c (x) = 1 outside of 0, then (|6.2p - (|6.4p imply 

Uu = Au,{x,t)e (E^Xfi) X (o,r), 
u (x, 0) = Ut (x, 0) = 0, X G IR^\ri, u\sj,= g (x, t) . 

Solving this initial boundary value problem in the domain {(x,t) G (]R'^\f2) x (0,r)} , we uniquely obtain 
the Neumann boundary condition p (x, t) for the function u, 

dnu\sT= P{x,t) ■ (6.5) 

CIPs arc quite complex problems. Hence, to handle them, one naturally needs to impose some simplifying 
assumptions. In this particular CIP our theory of the adaptivity is not working unless we replace the 
(5— function in (|6.3p by a smooth function, which approximates 5 {x — xq) in the distribution sense. Let 
>c G (0, 1) be a sufficiently small number. We replace S {x — xq) in (|6.3p with the function (x — xq) , 

(x - xo) = I ^-^"P ' " < [s^ix-xo)dx = l. (6.6) 

1 0, |x — xol > X, J 

We assume that k is so small that 



5>c (x — Xo) — in VL. 
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(6.7) 



Wc now introduce state and adjoint problems. Let € (0, 1) be a sufficiently small number. Consider the 
function G C°° [0,T] such that 

r [0,T-2C], 
zc{t) = { o,<e[r-C,T], (6.8) 

[ between and 1 for t e [0, T - 2C, T - C] • 

State Problem. Find the solution v {x, t) of the following initial boundary value problem 

c {x) Vtt ~ Aw — in Qt, 

w(a;,0) = i;t(a;,0) = 0, (6.9) 

Adjoint Problem. Find the solution A {x, t) of the following initial boundary value problem with the 
reversed time 

c {x) Xtt - AA = in Qt, 

X{x,T) = \t{x,T) = Q, (6.10) 

dnX \st H W id - (^1 • 

Here functions v e (Qt) and X £ (Qt) are weak solutions of problems (|6.9p and (j6.10p respectively. 
In fact, we need a higher smoothness of these functions, which we specify below. In (|6.9p and (|6.10p functions 
g and p are the ones from (|6.4p and (|6.5p respectively. Hence, to solve the adjoint problem, one should solve 
the state problem first. The function {t) is introduced to ensure the validity of compatibility conditions 
at {t = T} in (|6?T0l) . The Tikhonov functional for the above CIP is 

Eaic) = i J (vIst ^ g^x, t)YzQ (t) dadt + J (c - Cgiobfdx, (6.11) 

St n 

where the function Cgiob S C (fi) is the approximate solution obtained by our approximately globally con- 
vergent numerical method on the first stage of our two stage numerical procedure (section 1) and a is the 
small regularization parameter. 

State and adjoint problems are concerned only with the domain 51 rather than with the entire space . 
We define the space Z as 

3 

z = {/ : / e n i/i (11) ,c,. &L^{n),i = 1,2,3} , ll/IU = ll/llc(n) +E • 

i=l 

Clearly H <Z Z as a, set. To apply the theory of above sections, we express in subsection 6.2 the function 
c(x) via standard piecewise linear finite elements. Hence, we assume below that c € V, where 

r = {c e Z : c e (1 -a;,d + a;)}. (6.12) 

To find the Frcchet derivative of the functional Ea{c), wc need to find Frechet derivatives of functions 
solutions v,X of problems (|6.9p . (|6.10p . This, in turn requires a higher smoothness of functions p^g [111 I14j . 
Theorem 6.1 can be easily derived from a combination of Theorems 4.7.1, 4.7.2 and 4.8 of [11] as well as 
from Theorems 3.1, 3.2 of [14] . 

Theorem 6.1. Let ft C MP be a convex bounded domain with the boundary dfl G and such that 
there exists a function a £ (Tl) such that a \dn~ 0,d„a \dn= 1. Assume that there exist functions 
P {x, t) , $ (x, t) such that 

PeH^ (Qt) , $ e i?' (Qt) ; d^P \st= p (x, t) , |st= (*) 9 {x, t) , 
d(P{x,0) = dl<P{x,0) = 0,j = 1,2,3,4. 
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Then for every function c £ Y functions v, A G (Qt) i where v, A are solutions of state and adjoint 
problems \6.9\) . \6.1U\) . Also, for every c £ Y there exists Frechet derivative E'^{c) of the Tikhonov functional 
Ea-.Y ^Rin ifOI]) and 

T 

E'a{c) {x) = a (c - Cglob) {x) - j (utXt) {x, t) dt a (c - Cgiob) {x) +y{x) . (6.13) 



Functions E'^ic) [x) , y {x) £ C (fi) and there exists a constant D ~ D {Q, a, d, ui, z^) > such that 

Mc(n) < \\4l(n)^^PiDT) (||P||?,«(Q,) + ml^Qr)) ■ (6-14) 
The functional of the Frechet derivative E'^ (c) acts on any function b £ Z as 

E'^ic) ib)= J E'^ic) {x)b{x)dx. 
n 

6.2. Relaxation property for the functional Ea{c). In this section we use Tlieorems 5.2, 5.4 to 
derive tlie relaxation property for tire for the specific functional E^ (c) for our CIP. The first step is to define 
the operator F for our specific case. Set G := Y D H. We consider the set G as the subset of the space H 
with the same norm as the one in H. In particular, G = {c{x) £ H : c{x) £ [1 — w, d + w] for a; G Si} . Let 
H2 := L2 (St) ■ We define the operator F as 

F -.G ^ H2,F{c){x,t) ^ z^{t)[g{x,t)~v{x,t,c)], {x,t) & St, (6.15) 

where the function v := v{x,t,c) is the weak solution (??) of the state problem (|6.9p . g is the function 
in (|6.4p and {t) is the function defined in (|6.8p . For any function b E H consider the weak solution 
u (x, t, c, 6) e (Qt) of the following initial boundary value problem 

c (x) utt ^ Au-b (x) vtt, {x, t) £ Qt, 
u {x, 0) = ut {x, 0) = 0, u \st= 0. 

Theorem 6.2 can be easily derived from a combination of Theorems 4.7.2 and 4.10 of [TT]. 

Theorem 6.2. Let C 6e a convex bounded domain with the boundary dfl £ . Suppose that there 
exist functions a (x) , P (x, i) , $ [x, t) satisfying conditions of Theorem 6.1. Then the function u (x, t, c, b) £ 

{Qt) ■ Also, the operator F in H6.15\) has the Frechet derivative F' (c) (6) , 

F'(c) (6) = {t)u{x,t,c,b) |sr,Vce G,yb£H. 
Let B = B (i7, a, d, lo, zq) > be the constant of Theorem 6.1. Then 

\\F' (c)ll^ < exp {CT) (||P||ffa(Q,) + mHHQr)) , Vc G G. 
In addition, the operator F' (c) is Lipschitz continuous, 
\\F' (ci)-P'(c2)||^<exp(CT)(||P|| 

Following (|2.2p . we introduce the error of the level 6 in the data g{x,t) in (|6.4p . So, we assume that 

g{x,t) = g*{x,t) + gs{x,t); g* , gs £ L2 (St) ,\\gs\\L^^ST} - ^- (^'l^) 

where g*{x,t) is the exact data and the function gs{x,t) represents the error in these data. To make 
sure that the operator F is one-to-one, we need to refer to a uniqueness theorem for our CIP. However, 
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uniqueness results for multidimensional CIPs with single measurement data are currently known only under 
the assumption that at least one of initial eonditions does not equal zero in the entire domain f2, which is 
not our case. All these theorems were proven by the method, which was originated in 1981 in three papers 
[HllOlIig; also see, e.g. EH US |30l IH |3l|35] as well as sections 1.10, 1.11 of the book [11] and references 
cited there for some follow up publications of those authors about this method. This method is based on 
Carleman estimates. Although many other researchers have published about this method, we do not cite 
those works here, because the topic of uniqueness is not a focus of the current paper. We refer to surveys 
[551 H5] for more references. Lifting the above assumption is a long standing and well known open question, 
see [31] for a recent partial answer to this question. Nevertheless, because of applications, it makes sense 
to develop numerical methods for the above CIP, regardless on the absence of proper uniqueness theorems. 
Therefore, we introduce Assumption 6.1. 

Assumption 6.1. The operator F (c) defined in i6.15]) is one-to-one. 

Theorem 6.3 follows from Theorems 3.3, 6.1 and 6.2. Note that if a function c G is such that c G [1, d] , 
then by c e G. 

Theorem 6.3. Let D, C be a convex bounded domain with the boundary dfl G . Suppose that 
there exist functions a (x) , P (x, i) , $ (x, satisfying conditions of Theorem 6.1. Let Assumption 6.1 and 
condition \6.16]) hold. Let the function v ~ v {x, t, c) G (Qt) in ^6.11\) be the solution of the state problem 
h6.9fl for the function c € G. Assume that there exists the exact solution c* G G, c* (x) G [1, rf] of the equation 
F (c*) =0 for the case when in H6.16\) the function g is replaced with the function g* . Let in i6.16\) 

a^a{d)^ (5^^,^^ = const. G (0,1/4). 

Also, let in fi6.ll]) the function Cgiob G G and 

WCglob - C*|| < — . 



Then there exists a sufficiently small number 5^ = 5q \^il,d,uj,Zi^,a, \\P\\ h'>{Qt) ' ll^ll-fr^(QT) ''^j ^ such 
that Vga^^ (c*) C G,V(5 G (0,(5o) and the functional Ea (c) is strongly convex in Vg3^, (c*) with the strong 
convexity constant a/4. In other words, 

\\ci - C2II' < {E[, (ci) - E'^ (C2) , ci - C2) , Vci,C2 G Vs.. (c*) , (6.17) 

where (, ) is the scalar product in L2 {^) and the Frechet derivative E'^ is calculated via i6.13]) . Furthermore, 
there exists the unique regularized solution Cq,(5) , and c^^s^ G V^sn ^3 (x*) . In addition, the gradient method of 
the minimization of the functional Ea (c) , which starts at Cgiob, converges to Ca(s) - Furthermore, let ^ G (0, 1) 

be an arbitrary number. Then there exists a number 5i = 5i {vi,d,UJ, Zi^,a, \\P\\H<^(Qr^) : Il*5'lli/5(g5,) i/^iC^ ^ 
(0, (5o) such that for all 5 G (0, (5i) 



1(5) - C II < 



i \\cgiob - c*|| , if Cgiob 7^ c*, 

^, if Cgiob = C*. 



In other words, the regularized solution Ca[s) is more accurate than the solution obtained on the first stage 
of our two-stage numerical procedure. Furthermore, since E'^^^-^ (cq(5)) — 0, then {6.17^ implies that 



II 2 

^(-5)11 ^ ^ 



K(5) (c) ^ ,VcG (c*). 



Theorem 6.4 follows from Theorems 5.1 and 6.3 as well as from Theorem 4.11.3 of jllj . 
Theorem 6.4. Let conditions of Theorem 6.3 hold. Let \\c*\\ < B, where the constant B is given. Let 
M„ d H be the subspace obtained after n mesh refinements as described in section 4. Let h„ be the maximal 
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grid step size of the subspace M„. Let D ~ D {fl,a,d,uj^ Zi^) > be the constant of Theorem 6.1 and K be 
the constant in I4.8\j . There exists a constant N2 = N2 (^D,T, WPWn^'iQT) ' ll*^ll-ff^(QT)) that if 



5BN2K 



then there exists the unique minimizer c„ of the functional i6.11]) on the set G H M„ ■ Furthermore, Cn G 
Vg3^i (x*) n Mn and the following a posteriori error estimate holds 



II 2 
="(5)11 < 



K,s, (c.) . (6.18) 



The estimate (|6.18p is a posteriori because it is obtained after the function c„ is calculated. Theorem 
6.5 follows from Theorems 5.2, 5.3, 6.4, also see Theorem 4.11.4 in [TT] . 

Theorem 6.5 (relaxation). Assume that conditions of Theorem 6.4 hold. Let c„ G Vgsi^ (x*) D Mn be 
the unique minimizer of the Tikhonov functional fid. 11]) on the set GO A/„ (Theorem 6.4). Assume that the 
regularized solution c^^^^) 7^ c„, i.e. c^^^^ ^ A/„. Let rj e (0, 1) be an arbitrary number. Then one can choose 
the maximal grid size hn+i ~ hn+i (^B, N2, K,S, z^^, ij.,ri^ G (0, /i„] of the mesh refinement number {n + 1) 
so small that 

\\cn+l ~ Ca{S}\\ < V \\cn " Ca(5)\\ , (6.19) 

where the number N2 was defined in Theorem 6.4- Let ^ € (0, 1) be an arbitrary number. Then there exists 
a sufficiently small number 5a = 5o (^A, N2, K,S, Z(^,£_, ii,r]^ G (0,1) and a decreasing sequence of maximal 

grid step sizes {/ifc}^^^ = h^. [B, N2, K, 6, zq, ^, ii.rf) such that if 5 E (0, 5q) , then 

\\ck+i - c*\\ < r;^ - c„(,)|| + ( ^ W^^ioy , if c ^ c*, ^ ^ ^ _ 

" wil 1^ it Cgjgj = c , 

Theorem 6.6 follows from Theorems 5.4 and 6.5. 

Theorem 6.6. (relaxation for local mesh refinements). Assume that conditions of Theorem 6.5 hold. Let 
f2 = r2iUr22. Suppose that mesh refinements are performed only in the subdomain 0,2. Let /i'-^-' be the maximal 
grid step size in Vli. Then there exists a sufficiently small number 5q = (5o (^B, N2, K, z^, ij,,ri^ £ (0, 1) and a 

hkj jh-k — hk {B, N3, K, Z(^, fi, r/) of meshes in 02 such 

that if ||Vcq(5)||^ jj-j^j is so small that if 

ll^'^"('5)llL^(ni) ^^^^ - I 11'='^ ~ '="('5)11 l,...,n and 5 G (0, Jo) , 
then \6.20\l holds with the replacement of {hkYk=i V^/^ 1 

^ J k—l 

Here the number N3 = N3 (^D, T, ||P||ff6(Q^) , ll*!lif5(Q^)) > 0. 

7. Mesh Refinement Recommendations and the Adaptive Algorithm. 

7.1. Mesh Refinement Recommendations. Recommendations for mesh refinements are based on 
the theory of section 6. We now present some partly rigorous and partly heuristic considerations which 
lead to these recommendations. The latter means that both mesh refinement recommendations listed below 
should be verified numerically. We come back to the arguments presented in the paragraph above Theorem 
5.4. To simplify the presentation, assume, for example that 



Vcq(5) (x) = Vc* (a;) = for x eOi. 
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(7.1) 



A more general case when functions Ca(s) (x) ,c* (x) change slowly in fli can be considered similarly. Using 
(|4Jl) and (fTT]) . we obtain that {ca{s) ~ PkCa{S}) (x) = for a: e fii,Vfc > 1. Hence, by (jlJ]) 



||ca(5) -P„+lC„(5)||^^(fj) - \\Cai5} ^ ^n+lCa(5) 1 1 (^^ ) < ^ 1 1 1 1 (O2 ) ^^"+1' 

where hn+i is the maximal grid step size in fl2 after n + 1 mesh refinements. Hence, using the second equality 
((E7)) and we obtain 



I II 

|C„+1 -C„(5)|| < 



■^a(5)(Cn+l) ||^Ca(5)||^^^j^^^ /l„+l. (7.2) 



Given a function f E C (fl) , the main impact in the norm ||./||^2(f2) is provided by neighborhoods of those 
points X G where the function |/(a;)| achieves its maximal value. Hence, (j7.2p indicates that we should 
decrease the maximal grid step size /i„+i (i.e. refine mesh) in neighborhoods of those points x (z where 
the function (c„+i) (.t)| achieves its maximal values, where the function E'^^^^ (cn+i) (a^) € C (fJ) is given 
by formula (j6.13p . Although after n mesh refinements we know only the function c„ € M„ rather than the 
function c„+i € M„+i, still, since functions c„ and c„+i are sufficiently close to each other, we should likely 
refine mesh in neighborhoods of those points x G where the function \E'^ (c„) (a;)| achieves its maximal 
values. These considerations lead to two mesh refinement recommendations below. 

The First Mesh Refinement Recommendation. Let f3i e (0, 1) be the tolerance number. Refine 
the mesh in such subdomains of 0,2 where 

\K{Cn){x)\>Pl'^^\E'ACn){x)\. (7.3) 

To figure out the second mesh refinement recommendation, we note that by (|6.13p and (j6.14p 

K(S) (cn) {x) < a {\\cn\\c(n) + kaiob\\c(n)) + l|c„|lc(a) exp(i:'T) (!|P|Ih6(q^) + IMI^q^)) ■ 

Since a is small, then the second term in the right hand side of this estimate dominates. Next, since we have 
decided to refine the mesh in neighborhoods of those points, which deliver maximal values for the function 



a(5) ('^") (^) 



then we obtain the following mesh refinement recommendation. 

Second Mesh Refinement Recommendation. Let € (Oj 1) the tolerance number. Refine the 
mesh in such subdomains of O2 where 

Cn [x) > P2 maxcn (a;) , (7.4) 

In fact, these two mesh refinement recommendations do not guarantee of course that the minimizer 
obtained on the corresponding finer mesh would be indeed more accurate than the one obtained on the 
coarser mesh. This is because right hand sides of formulas (|7.3p and (|7.4p are indicators only. Thus, numerical 
verifications are necessary. As to tolerance numbers Pi and /32, they should be chosen numerically. Indeed, 
if we would choose /3i,/32 ~ 1; then we would refine the mesh in too narrow regions. On the other hand, 
if we would choose fii,P2 ~ 0, then we would refine the mesh in almost the entire subdomain f22, which is 
inefficient. 

7.2. The adaptive algorithm. Since this algorithm was described in detail in a number of publica- 
tions, see, e.g. [HI [14], we outline it only briefly here. Recall that the adaptivity is used on the second stage 
of our two-stage numerical procedure (section 1). On the first stage the approximately globally convergent 
algorithm is applied. It was proven, within the framework of the so-called Second Approximate Mathemati- 
cal Model, that this algorithm delivers a good approximation for the exact solution c* [x) of the above CIP, 
see Theorem 2.9.4 in [TT] as well as Theorem 5.1 in [T7]. We start the adaptivity on the same mesh on which 
the algorithm of the first stage has worked. In our experience, this mesh does not provide an improvement 
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(a) Gfdm 



(b) G — GpEhi U Gfdm 



(c) Gfem ^ 



Fig. 8.1: The hybrid mesh (b) is a combinations of a structured mesh (a), where FDM is applied, and a mesh (c), 
where we use FEM, with a thin overlapping of structured elements. The solution of the inverse problem is computed 
in the square and c{x) = 1 for x G G\^. 



of the image. On each mesh we find an approximate solution of the equation E'^ (c) = 0. Hence, by (|6.13p 
we find an approximate solution of the following equation on each mesh 



For each newly refined mesh we first linearly interpolate the function Cgiob [x) on it. Since this function was 
initially computed as a linear combination of finite elements forming the initial mesh and since all our finite 
elements are piecewise linear functions, then subsequent linear interpolations on finer meshes do not change 
the function Cgiob (x). On each mesh we iteratively update approximations cJJ of the function 0^(5). To do 
this, we use the quasi- Newton method with the classic BFGS update formula with the limited storage |42j . 
Denote 



where functions Vh (x, c^) , \h (a;, t, d^) ^''"C FEM solutions of state and adjoint problems ()6.9p . (|6.10|) with 
c := c^. We stop computing cJJ if either ||<p"||L2(n) < 10^^ or norms 1 1</'"| liain) are stabihzcd. Of course, 
only discrete norms are considered here. 

For a given mesh obtained after n mesh refinements, let c„ be the last computed function on which 
we have stopped. Next, we compute the function \E'^ (c„) {x)\ using (|6.13p . where v := Vh {x,t,Cn) , A := 
Xii {x,t, Cn) ■ If we use both above mesh refinement recommendations, then we refine the mesh in neigh- 
borhoods of all grid points satisfying (|7.3p and (|7.4p . In some studies, however, we use only the first 
recommendation. In this case we refine the mesh in neighborhoods of all grid points satisfying only (j7.3p . 

8. Numerical Studies. We present here three numerical examples of the performance of our two-stage 
numerical procedure: one for computationally simulated and two for experimental data. More numerical 
tests of the adaptivity technique can be found in[ll|6l[7l[8l[9l[l0l[ni[l3l[l4l[l5l[16]. In Test 1 we have used 
only the First Mesh Refinement Recommendation, and in Tests 2,3 we have used both recommendations. 
Since the numerical method of the first stage of our procedure is not a focus of this paper, and since it was 
described earlier in, e.g. [TTJ [T31 [HI [HI [311 [371 [35], we do not describe it here. 



T 
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8.1. Computationally simulated data. Test 1. Wc conducted computational simulations in two 
dimensions. Since it is impossible to computationally solve equation (|6.2I) in the entire space R^, we have 
conducted data simulations in the rectangle G = [—4, 4] x [—5, 5] . To simulate the boundary data g {x, t), we 
have solved the forward problem by the hybrid FEM/FDM method [5] using the software package WavES 
[47] . To do this, we split the domain G in two subdomains G = Gfem U Gfdm, see Figure [01 Here 
Gfem '■— ^ = [—3,3] X [—3,3] and Gfdm = G\Gfem ■ The coeffieicnt c{x) is unknown in the domain 
ft C G and is defined as 

{1 in Gfdm, 

l + b{x) in Gfem, (8.1) 
4 in small squares of Figure 18. H 

where the function b{x) £ D,fem is defined as 



b{x) 



for(xi,a;2) G Qfem : -2.875 < xi < 0,-2.875 < < 0, 



0.5sin2(f|^)sin2(f|^) otherwise. 



The spatial mesh consists of triangles in Gfem and of squares in Gfdm with the grid step size h = 0.125 
both in overlapping regions and in Gfdm- There is no reason to refine mesh in Gfdm since c{x) = 1 in 
Gfdm- Let dGi and 96*2 be, respectively, top and bottom sides of the rectangle G and dGs be the union of 
vertical sides of G. Wc use first order absorbing boundary conditions on dGi U 96*2 [24] and zero Neumann 
boundary condition on dG^. 

Let s be the upper value of the Laplace transform of the solution of our forward problem. We use this 
transform on the first stage of our two-stage numerical procedure. It was found that for the above domain 
n the optimal value is s = 7.45. Consider the function / (t) , 



fit) 



0.1 [sin {st - 7r/2) + 1] , t G [0, ti] , ti = 27r/s, 
0,te {ti,T],T= 17.8ti. 



The forward problem for data simulations is 

c(x)m„-Au = 0, inGx(0,r), 

0) = 0, ut{x, 0) = 0, in G, 

^""laci ^ ^9tu, on dGi x (ti,T), ^^''^^ 

a„u|gg^ = 0, on dG3 X (0,T). 

The solution of this problem gives us the function g{x,t) = u\st - Next, the coefficient c{x) is "forgotten" 
and we apply the two-stage numerical procedure to reconstruct it from the function g {x,t) . To have noisy 
data, we have added the random noise to the function g (x, t) as 

9i,3 = 9 {x',t^) [1 + 0.02aj (gmax - ffmin)] • (8.3) 

Here € dil and P S [0,T] are mesh points on dft and [0,r] respectively, g^i^ and 5max are minimal and 
maximal values of the function g and aj € [—1,1] is the random variable. The "inverse crime" was not 
committed here since we have introduced the noise in the data and because the grids in both stages of our 
two-stage numerical procedure were different from the one which was used to solve the problem (|8.2[) . 

1. The approximately globally convergent stage. Since we focus on the adaptivity in this paper, we do 
not describe this algorithm here and refer to section 2.6.1 of [TT] instead. Figure [5?2] displays the result of 
this stage. 
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X-Dl splaf«n 
H zaazi 

^ 2.441 ( 

■ 2.201 E 
1.951; 

a) exact coefBcient c{x) b) CoefBcient Cgiot reconstructed on the first stage 

Fig. 8.2: a) Spatial distribution of the exact coefRcient c{x). b) Result of the performance of the approximately 
globally convergent algorithm (first stage). The spatial distribution of the computed coefficient Cgiob displayed. Here 
maxCgioi, (x) = 3.2, whereas maxc(x) — 4. Hence, we have 20% error in imaging of the maximal value of the function 
c (x) . The slowly changing part of the function c (x), i.e. the second raw in the above definition of the function b (x) , 
is not imaged. Comparison with Figure a) shows that while the location of the right inclusion is imaged correctly, 
the left one still needs to be moved upwards. This is done on the second stage of our two-stage numerical procedure, 
i.e. on the adaptivity stage. On this stage we take the function Cgiob (x) as the starting point for the minimization of 
the Tikhonov functional (|6.11[l . The second stage refines the image of the first. 



2. The adaptivity stage. Since we have observed that u {x, T) « 0, we have not used the function 
(t) in our computations. In this test we take the noise level 2% in (|8.3p and the regularization parameter 
a ~ 0.02 in (j6.1ip . We now comment on the stopping criterion for mesh refinements, which we use in 
numerical studies of the adaptivity technique in this paper. Let c„ is the coefficient c{x) calculated after 
n mesh refinements. In Theorems 5.2-5.4, 6.5, 6.6 the relaxation parameter ry is independent on the mesh 
refinement number n. In practice, however, one should expect such dependence 77 := 77„. In this case the 
parameter ry of those theorems is 77 = max(r7„). Then because of the relaxation property of Theorems 6.5, 
6.6 as well as because of Remark 5.1, it is anticipated that numbers 77„ decrease with the grow of n until 
the regularized solution Ca(s) is approximately reached. However, nothing can be guaranteed about numbers 
?/„ as soon as the regularized solution is reached. Hence, in our computations of the adaptivity method we 
stopped mesh refinement process at such n := no that rj^^ > ?7„q_i. If rj^^ « Vno-i^ then we took the final 

solution C final ■= Cno- 

Figure [S751 -e) . f) represents the images obtained after 4 and 5 mesh refinements, respectivelly, as well as 
adaptive locally refined meshes are presented on l8.3l -a)-d). Comparing with Figure [Pi c), one can observe 
that locations of both inclusions are imaged accurately. Recall that in each inclusion of Figure lS.U c) c {x) = 4, 
see definition for c{x) in (|8.ip shown also on Figure [ST^ a) . Therefore, maximal values of the function c{x) 
on Figures [5T51 -e).f) are also accurately imaged: the error docs not exceed 3.5%. 

Figure [8T4l displays the graph of the dependence of the norm ||c„ — Calj^^^^fj-) from the mesh refinement 
number n. By (|6.19p and (|6.20p these norms should decay. Since we do not exactly know what the regularized 
solution Ca is, we have taken Cq := C4 on Figure [8^ a). On Figure IS^ b) wc have superimposed those graphs 
for Ca '■— C4 and Ca '■= C5. One can observe that norms ||c„ — Ca\\ decay in the case when Ca is taken on 
the 4-th refined mesh. At the same time we also observe, that the relaxation property (|6.20p is not fuUfiUed 
when we take Ca on the 5-th refined mesh since 773 > 772, see l8.4l -b). Thus, we take the final reconstruction 
Ca ■= C4, the function obtained after four (4) mesh refinements. 

Remark 8.1. It is well known that imaging of locations of small inclusions and maximal values of the 
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a) 4776 elements 



b) 5272 elements 



c) 6174 elements 




d) 7682 elements e) C4(x), max 04(0;) = 3.9 f) C5(x), max 05(2::) = 3.87 



Fig. 8.3: Adaptively refined meshes (a)-(d) and finally reconstructed images (e) and (f) on 4-th and 5-th adaptively 
refined meshes, respectively. On e) maxc4 = 3.9 and on f) maxcs — 3.87. Reconstructed function on e) is obtained 
on the mesh presented on d). The mesh for the function on f) is not shown. Locations of both squares of Figure 
18.31 -a) as well as maximal values of the computed funtion Cgiot (x) in them are imaged accurately. 




a) b) 

Fig. 8.4: a) Computed relaxation property ||c„-|-i — Cq||i,, < 77„||c„ — CqUlj for ^-he noise level 2% in (|8.3|l and the 
regularization parameter a — 0.02 in (|6.1ip . Here, < 77^ < 1 is the small relaxation parameter obtained after n 
mesh refinements. Here, we take Ca on the 4-th refined mesh shown on the Figure I8.31 d') . b) Comparison of the 
relaxation property ||c„+i — Cc«|!l2 < 'ynllcii — Ca||L2 when we take different functions Ca'. on the 4-th or on the 5-th 
refined mesh. 
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Source 




Fig. 8.5: Schematic diagram of data collection. Original source: M. V. Klibanov, M. A. Fiddy, L. Beilina, N. Pan- 
tong and J. Schenk, Picosecond scale experimental verification of a globally convergent numerical method for a coeffi- 
cient inverse problem, Inverse Problems, 26, 045003, doi:10. 1088/0266-5611/26/4/045003, 2010. ©lOP Publishing. 
Reprinted with permission. 



function c {x) in them is of the primary interest in apphcations and it is more interesting than imaging of 
slowly changing parts. Indeed, small inclusions can be explosives |371 138] . tumors, etc.. 

Remark 8.2. The above stopping criterion for mesh refinements shows that relaxation Theorems 6.5, 
6.6 are quite useful for computations. 

8.2. Experimental data. Experimental studies were described in detail in [TS1[32] as well as in Chapter 
5 of . Hence, we omit many details here. We point out that the main difficulty was a huge misfit between 
computationally simulated and experimental data. The latter was the case even for the free space data: the 
analytic solution predicted by Maxwell equations was radically different from the experimentally measured 
curves. This can be explained by unknown nonlinear processes in both transmitters and detectors. The same 
was observed for the backscattcring data collected in the field, see [35] and section 6.9 of [TT]. To handle this 
misfit, a new data pre-processing procedure was applied. This procedure has immersed experimental data 
in computationally simulated ones, see Figures 4 in |38j and Figures 5.3 in Naturally, this procedure 
has introduced a significant modeling noise in already noisy data. Nevertheless, computational results were 
very accurate ones, which speaks well for the robustness of our reconstruction method. The first stage of our 
two-stage numerical procedure was working with blind data (unlike the second stage). Therefore, results of 
at least the first stage were unbiased. 

The data collection scheme is displayed on Figure [875] A single source of electric wave field emits pulse 
for only one component of the electric field, two other components were not emitted. The prism is our 
computational domain il. The outcome time resolved signal was measured at many detectors located on the 
bottom side of the prism. The same component of the electric field was measured as the one emitted. Since 
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Case number 


Computed n 


Directly measured n 


Computational error 


1 (Cube 1) 


1.97 


2.07 


5% 


2 (Cube 1) 


2 


2.07 


3.4% 


3 (Cube 1) 


2.16 


2.07 


4.3% 


4 (Cube 1) 


2.19 


2.07 


5.8% 


5 (Cube 2) 


1.73 


1.71 


1.2% 


6 (Cube 2) 


1.79 


1.71 


4.7% 



Table 8.1: Blindly computed and directly measured refractive indices n by the first stage of our two-stage numerical 
procedure. The error in direct measurements was 11% for cases 1-4 (Cube 1) and 3.5% for cases 5,6 (Cube 2). 



we have not measured that signal at the rest difl of 9f2, we have prescribed to diQ the same boundary 
conditions as ones for the uniform medium with the dielectric constant = 1. The prism Q is filled with 
a dielectric material with the dielectric constant Sr ~ 1, i.e. almost the same as in the air. We point out, 
however, that when using the first stage of our two-stage numerical procedure, we did not use any knowledge 
of the dielectric constant of this prism. We have only used the fact that Sr = I outside of this prism, see 

(EH). 

Wc have placed one dielectric inclusion inside of this prism. Inclusions were two wooden cubes, which we 
call below "Cube 1" and "Cube 2" . Sizes of their sides were 4 cm for Cube 1 and 6 cm for Cube 2. Note that 
only refractive indices n = rather than dielectric constants can be measured directly in experiments. 
The goal of the first stage was to reconstruct the refractive index of the inclusion and its location. The 
goal of the second stage was to reconstruct all three components of inclusions: refractive indices, shapes 
and locations. Since only one component of the electric field was measured, we have modeled the wave 
propagation process via the problem (|8.2p with Sr {x) := c (x), where the domain G C M'^ was a prism, which 
was bigger than the prism ft, see (5.8) and section 5.4 in [11] for this domain. The function / (t) in (|8.2p was 

r sin {ujt) ,te{0, 2n/uj) , 
^^^"\ 0,t>27r/w, 

where a; = 14 for Cube 1 and a; = 7 for Cube 2 (see page 329 of [TT] and page 26 of [12] for w). It was 
only later, after the first author has conducted numerical simulations for solving the Maxwell equations |18| . 
when we have realized that the choice of modeling by one PDE only was well justified. In our experiments. 
Cubes 1 and 2 were placed total in six different positions. 

Table 18.11 summarizes results of blind study of the first stage of our two-stage numerical procedure. 
Because of the blind test requirement, direct measurements of refractive indices were performed by the 
conventional so-called "waveguide method" [H] only after computations of the first stage were done. Next, 
computational results were compared with measured ones. One can see that we had only a few percent 
difference with a posteriori directly measured refractive indices of both cubes. Furthermore, in five out of 
six cases this error was even less than the error in direct measurements. 

We now focus on the results which we have obtained on the second stage of our two-stage numerical 
procedure when applying the adaptivity. 

Test 2. The two stage numerical procedure for Case 1 of Table [8711 Figure [5751 -a) displays 
the result of the first stage of the two-stage numerical procedure. One can see that although the refractive 
index n ~ 1.97 and location of the inclusion are accurately calculated, the shape is inaccurate. The image 
of Figure [8^ a) was taken as the starting point for the adaptivity technique for refinement. The result of 
the second stage is presented on Figure IH^ b*). One can see that all three components of the inclusion are 
accurately reconstructed. In addition, the values of the function Sr (x) := c (x) ~ 1 outside of the imaged 
inclusion are also accurately computed. 

Test 3. The two stage numerical procedure for Case 6 of Table [8711 Figures [5771 -a) and 18.71 -b) 
display computational results for first and second stages, respectively. The rest of comments are the same 
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i 




a) max Ugiob = max ^er,giob 



max^c^ 



1.97 



b) maxn 



max yfe^ ~ max ^/c ~ 2.05 



Fig. 8.6: Case 1 of Table ISTTl was tested by the two-stage numerical procedure, a) The computational result of the 
first stage. Both location of the inclusion and refractive index Ugiob = 1-97 are accurately reconstructed. However, 
the shape of the inclusion is not reconstructed accurately, b) The computational result of the second (refinement) 
stage. All three components of the inclusion are very accurately reconstructed: refractive index, location and shape. 
Also, values of the function Sr (x) :— c (x) — 1 outside of the imaged inclusion are computed very accurately. Original 
source: L. Beilina and M.V. Klibanov, Reconstruction of dielectrics from experimental data via a hybrid globally 
convergent/adaptive inverse algorithm. Inverse Problems, 26, 125009, doi:10. 1088/0266-5611/26/12/125009, 2010. 
©lOP Publishing. Reprinted with permission. 



as ones for Test 2. Note that the shape is now reconstructed better than in Test 2. This can be heuristically 
explained as follows. The wavelength of our electromagnetic wave was /i = 3 cm. Thus, the size of the side 
of Cube 1 is 4 cm=1.33^. One the other hand, the size of the side of Cube 2 is 6 cm=2^, which is larger. 
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a) maxngiob = max y/Sr^giob ~ max y/Cgiob ~ 1-79 b) maxn — maxy^ ~ maxy^ = 1.73 

Fig. 8.7: Case 6 of Table [871] was tested for the two-stage numerical procedure, a) The computational result of the 
first stage. Both location of the inclusion and refractive index Hgiob = 1-79 are accurately reconstructed. However, 
the shape of the inclusion is not reconstructed accurately, b) The computational result of the second (refinement) 
stage. All three components of the inclusion are accurately reconstructed: location, refractive index, and shape. In 
addition, values of the function Sr (x) := c (a;) = 1 outside of the imaged inclusion are computed accurately. Original 
source: L. Beilina and M.V. Klibanov, Reconstruction of dielectrics from experimental data via a hybrid globally 
convergent/adaptive inverse algorithm. Inverse Problems, 26, 125009, doi:10. 1088/0266-5611/26/12/125009, 2010. 
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step 1 5 

Contour Fill of Scalar result. 



